A landscape of mouse mitochondrial small non-coding RNAs

Small non-coding RNAs (ncRNAs), particularly miRNAs, play key roles in a plethora of biological processes both in health and disease. Although largely operative in the cytoplasm, emerging data indicate their shuttling in different subcellular compartments. Given the central role of mitochondria in cellular homeostasis, here we systematically profiled their small ncRNAs content across mouse tissues that largely rely on mitochondria functioning. The ubiquitous presence of piRNAs in mitochondria (mitopiRNA) of somatic tissues is reported for the first time, supporting the idea of a strong and general connection between mitochondria biology and piRNA pathways. Then, we found groups of tissue-shared and tissue-specific mitochondrial miRNAs (mitomiRs), potentially related to the “basic” or “cell context dependent” biology of mitochondria. Overall, this large data platform will be useful to deepen the knowledge about small ncRNAs processing and their governed regulatory networks contributing to mitochondria functions.


Introduction
The recent advances in high-throughput sequencing technologies shed light on a plethora of RNA species with limited coding potential, collectively called non-coding RNAs (ncRNA) [1].Besides ribosomal RNA, ncRNA molecules are generally classified according to their size threshold in long ncRNA (lncRNAs), longer than 200 nt up to several kilobases, and small or short ncRNAs, from few to 200 nt.lncRNA are involved in the regulation of gene expression at transcriptional and post-transcriptional levels by interacting with DNA, proteins and other RNA molecules.Short ncRNAs include the well known tRNA (transfer RNA), involved in the translation of mRNA; snRNA (small nuclear RNA), involved in splicing; snoRNA (small nucleolar RNA), implicated in ribosomal RNA modification; Piwi-interacting RNA (piRNA) and microRNAs (miRNA) [2].piRNAs, 24-35 nucleotide long, 2'-O-methyl modified at 3'-end, associate with PIWI proteins; they were discovered in germlines in 2001, representing the most abundant small RNAs in mammalian testes and are essential for fertility; they are classically associated with "genome defense" against transposable elements and viruses [3][4][5]; more recently, they have been involved in epigenetic modification and gene expression regulation, probably at lower extent also in somatic tissues [6].piRNA biogenesis and functioning are still not completely known, probably due to the occurrence of different production pathways, the wide diversity of functions (in the nucleus and cytosol) and poorly conserved pathways specifically associated to features of germ cell development [6,7].
miRNAs, the most studied group of small ncRNAs, with a size of about 20 nt, negatively regulate gene expression by binding to target transcripts, thus affecting their translation and/or stability [8].In mammals, up to 50% of all protein-coding transcripts are predicted to be targeted by miRNAs; in addition, one miRNA can bind different mRNAs and one mRNA can targeted by different miRNAs, giving rise to complex regulatory networks that play key roles in almost all physiological pathways, but also in the pathogenesis of several diseases, especially cancer [9][10][11].miRNAs biogenesis is compartmentalized between nucleus and cytoplasm, but they are largely operative in the cytoplasm.However, emergent data indicate that mature miRNAs can localize in the nucleus and in subcellular compartments in the cytoplasm, unveiling possible new and specialized functions related to specific cellular compartments [12].In this regard, some studies have detected microRNAs in mitochondria isolated from rat, mouse and human cells, collectively called mitomiRs (mitochondria-localized microRNA) [13,14].In fact, in 2006 a few miRNAs were found in rat's liver mitochondria; although initially thought as cytosolic contamination, few years later other studies performed by microarray profiling and more recently by RNAseq confirmed the localization of miRNAs in mitochondria [15][16][17][18][19][20][21][22].
Mitochondria are the powerhouses in the cell and essential organelles for the integration of several key metabolic processes such as respiration, fatty and amino acid metabolism [23].They can be defined as central regulators of cellular homeostasis, and dysregulation of any activity has pathological consequences, like cancer, inflammation, neurodegeneration [24][25][26][27].Mitochondria contain their own circular genome, about 16 kbp in length, with 37 genes producing 2 rRNAs, 22 tRNAs and 13 proteins of the oxidative phosphorylation (OXPHOS) system [28].Mitochondrial homeostasis is dependent on the expression of its own genes, quite limited in number, but largely on the import of nuclear-encoded proteins from the cytoplasm [13].It is becoming increasingly clear that nuclear genome contributes to mitochondria activities not only through proteins (about 1500 nuclear-encoded proteins), but also through a number of miRNAs that can affect transcripts encoding mitochondrial-localized and mitochondrialrelated proteins.This crosstalk between nucleus and mitochondria is essential for regulating mitochondria dynamics (biogenesis, fusion, fission, mitophagy), metabolic pathways and mitochondria-mediated apoptosis, whose dysregulations have pathological consequences.Indeed, some miRNAs can be translocated to mammalian mitochondria.The proposed import mechanisms include a detachment of AGO2 and miRNA from RISC or just GW182 due to AGO2 phosphorylation or some other signal activation promoting their translocation together or separately; this process could be stimulated by P-bodies; translocation across mitochondrial membranes is unknown but suggested to be promoted by PNPase and occurs within TOM/TIM complexes, or relies on VDAC at the outer mitochondrial membrane [13,29,30].
Different gaps in the knowledge emerge from literature related to ncRNA contribution to mitochondria biology.As an example, mitomiR function still remains controversial, especially due to debated presence of RISC components required for their canonical mechanism of silencing activity [12,31].Moreover, a deep sequencing analyses of small RNA content in mitochondria, probably a foundation for a systematic understanding, is available only for HEK293 and HeLa human cells, and human myotubes [20,22], and regarding to mouse model only for primordial male germ cells, spermatogonia, spermatozoa, oocytes and zygotes [21].To gain insights into possible new and specialized function of small RNAs and in particular miRNAs related to mitochondrial biology, this work explores the content of small RNAs population by RNAseq of mitochondria isolated from different mouse tissues, also in order to identify a possible shared pool and distinctive mitomiR signatures potentially involved in basic or cell-context dependent mitochondria functions.

Animals and ethics statement
Three wild-type male CD1 mice, *3 months old, were used in this study.They were maintained in a temperature-controlled room at 22 ± 2˚C under a 12 h dark/light cycle with a standard pellet diet and unlimited access to water.
Mice were sacrificed under anesthesia by cervical dislocation.In detail, animals were placed in a plexiglas chamber with 4% isoflurane (Iso-Vet, Piramal Healthcare, United Kingdom) for 5 min and were sacrificed when fully sedated, as measured by a lack of heartbeat and active paw reflex.Tissues were rapidly removed, rinsed in PBS (pH 7.6) to remove blood contaminants and used for experimental procedures.Every effort was made to minimize animal pain and suffering.
Experiments were approved by the Italian Ministry of Education and the Italian Ministry of Health (authorization n˚48/2022-PR issued on 28.01.2022).Procedure involving animal care were carried out in accordance with the National Research Council's publication Guide for Care and Use of Laboratory Animals (National of Institutes of Health Guide).

Mitochondria purification
Mitochondria were purified from freshly excised tissues (liver, gastrocnemius skeletal muscle, testis, white adipose tissue) using Qproteome Mitochondrial Isolation kit (Qiagen) according to manufacturer's instructions with minor modifications and performing all the steps at 4˚C.Briefly, tissues were homogenized in Lysis Buffer supplemented with Protease Inhibitor solution by Potter homogenizer, incubated on an end-over-end shaker for 10 min and centrifuged at 1000 x g for 10 min.The pellet was resuspended in ice-cold Disruption Buffer supplemented with Protease Inhibitor, subjected to Potter homogenizer for cell disruption, followed by centrifugation at 1000 x g for 10 min.The supernatant was centrifuged at 6000 x g for 10 min to pellet down mitochondria.The mitochondria were further purified by resuspending the pellet in Purification buffer and carefully pipetting onto layers of Purification Buffer and Disruption Buffer.The solution was centrifuged at 14000 x g for 15 min.The mitochondrial ring at the interface of Purification Buffer/Disruption Buffer was collected and washed in Mitochondria Storage Buffer; the final mitochondrial pellet was resuspended in 20 μl Storage buffer.
In order to remove RNA of nuclear origin present in the cytosol and eventually adsorbed on the outer mitochondrial membrane, mitochondria were treated with RNase A solution, as already described [20,22].In brief, the 20 μl purified mitochondria from the above described procedure were mixed with 300 μl PBS containing RNase A (100 μg/ml) and incubated at 37˚C for 1 h.Rnase A treated mithocondria were pelleted down at 8000 rpm for 10 min, washed twice in PBS and finally resuspended in Qiazol buffer (Qiagen) for RNA extraction.

RNA extraction, control of mitochondria purity and real-time PCR analyses
Mitochondria and cytosolic fractions obtained by the procedure described above according to Qproteome Mitochondrial Isolation kit (Qiagen) protocol were used for total RNA purification by miRNeasy Mini Kit (Qiagen) according to the manufacturer's protocol.The yield and quality of RNA were evaluated by NanoDrop One spectrophotometer (NanoDrop), Qubit 4 Fluorometer (Invitrogen) and by TapeStation 4200 (Agilent Technologies).
The mitochondrial RNA enrichment and eventual cytosolic RNA contamination were checked by RT-qPCR.Mitochondria and cytosolic RNAs were retrotranscribed by SensiFAST cDNA Synthesis kit (Bioline).Then standard SYBR Green Real-time qPCR assays were performed for amplifying β-actin transcript as nucleus-derived RNA control and Cytochrome c as mitochondrion-derived control with the following primers: β-actin, 5'-CAACGGCTCCG GCAT -3' and 5'-CTCTTGCTCTGGGCC -3'; Cytochrome b, 5 0 -TACCTGCCCCAT CCAACATT-3 0 and 5 0 -TAAGCCTCGTCCGACATGAA-3 0 .For the further experimental steps, we selected the RNA mitochondrial preparations giving the following qPCR parameters on triplicate analyses: Ct>35 for β-actin, which was fully detectable in the parallel cytosolic sample, indicating the absence of cytosolic contaminant RNA; Ct<28 for cytochrome b, indicating a good enrichment of mitochondrial RNA.

RNA sequencing and data analyses
Indexed libraries were prepared from 50 ng of each purified RNA sample with NEXSmallRNA Seq v3.0 (PerkinElmer) according to the manufacturer's instructions.Libraries were quantified using the TapeStation 4200 (Agilent Technologies) and Qubit 4 Fluorometer (Invitrogen).Libraries were pooled in equimolar amounts with final concentration of 2 nM.The pooled samples were subject to cluster generation and sequencing using an Illumina NextSeq550Dx System (Illumina Inc.) in a 1x75 single-end format with an average depth of the sequences of 20 Million reads/sample.The raw sequence files generated (.fastq files) underwent quality control analysis using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/).sRNAbench by sRNAtoolBox was used to identify the non-coding RNAs signature in all sequenced samples, thus removing adapter sequences, based on reference kit instructions, and the low quality reads [32].In particular, the following steps were performed: 4 nt barcode from 5' end of reads were removed; the adapter sequence "TGGAATTCTCGGGTGCCAAGGG" were detected; minimum length of adapter sequence aligned considered as 10 and only 1 mismatch allowed between adapter and reference sequences; 4 nt random adapter from 3' end of adapter trimmed reads were removed.
Finally, reads were also aligned to mouse mitochondrial genome (NC_005089.1 from NCBI) by blastn command line applications with defaults parameters, finding reads with at least one reported alignment in the following percentages: 4.87% for gastrocnemius muscle, 6.17% for liver, 4.00% for testis and 3,47% for WAT.e-value�0.01was considered statistically significant.Finally, all piRNA sequences of piRbase (http://bigdata.ibp.ac.cn/piRBase/) and all miRNA sequences included in miRbase (https://www.mirbase.org/)were aligned to mitochondrial genome with parameters described above.

mitomiR targetome analyses
The list of 126 tissue-shared mitomiRs and the lists of tissue-specific mitomiRs were used to launch a search on miRWalk platform (http://mirwalk.umm.uni-heidelberg.de/) to retrieve all their target transcripts.In detail, by selecting "miRTarbase", or "TargetScan" and "miRDB", and "5 0 UTR", or "3 0 UTR" or "coding region" and selecting 0.5 as score for miRNA-mRNA pairings, 6 different target lists for tissue-shared and the 4 tissue-specific mitomiRs subset were retrieved containing all experimentally validated targets (by miRTarbase tool), or those predicted by both TargetScan and miRDB tools.By exporting the lists in Excel, a list of nonredundant transcripts representing the whole targetome of tissue-shared and 4 lists of tissuespecific mitomiRs were produced.On those lists, enrichment of biological pathways supplied by Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway annotation was performed by Database for Annotation, Visualization and Integrated Discovery (DAVID, https://david.ncifcrf.gov)which facilitates the transition from data collection to biological analysis by systematically extracting biological meaning and gene functional groups from large gene lists.Considering this type of analysis could increase the false-positive rate, to control the false-positive rate in the results, a multiple test correction of enrichment p values was performed on the functional annotation categories by selecting the Benjamini-Hochberg procedure.Biological pathways with p values < 0.05 were considered statistically significant [34].
Targetomes of tissue-shared or-specific mitomiRs were also analyzed by MitoCarta3.0 to extract enriched biological pathways specifically related to mithocondria biology.In brief, targetome lists were matched with 1140 mouse genes encoding proteins with strong support of mitochondrial localization; the merging target transcript lists, plus those mitochondrial transcripts reported in PubMed as experimentally validated target of our detected mitomiRs, were analyzed by MitoCarta3.0 for the assignment to 149 hierarchical 'MitoPathways' spanning seven broad functional categories relevant to mitochondria (http://www.broadinstitute.org/mitocarta) [35].

Rationale of the study
In order to have a large platform of data regarding the population of small RNAs in mitochondria, arguing their potential contribution to mitochondrial activities, novel mechanisms of regulation and crosstalk with other cellular components, we decided to select (Fig 1 ): i. liver, since it can be seen as the metabolic hub of the body, in analogy to mitochondria being the metabolic hub of the cell; in fact, the liver is extremely rich in mitochondria, with each hepatocyte containing 1000-2000 mitochondria [36]; ii. gastrocnemius skeletal muscle, whose health is dependent on the optimal function of its mitochondria [37]; due to the high metabolic requirements for locomotion, mitochondrial ATP production is very important, with typically 3-8% of the skeletal muscle volume being mitochondria [38]; iii. testis, since spermatogenesis relies on mitochondrial dynamics and functions, whose disruption causes male infertility [39]; iv.white adipose tissue (WAT), where mitochondrial activity is essential for the choice of carbon source (use of carbohydrate through glycolysis versus use of lipids for oxidative phosphorylation) and impacts on tissue and systemic energy homeostasis; the important role of mitochondria in the regulation of white adipose tissue remodeling and energy balance is increasingly appreciated [40].
Mitochondria were purified as described in Material and Methods section, and smallRNA sequencing experiments were performed on extracted RNAs (  ; their finding into mitochondria of testis is consistent with bioinformatic analyses on smallRNAseq libraries from primordial germ cells and spermatozoa of male mice [21,41].Surprisingly, here piRNAs were also detected in the mitochondria of the other tissues, representing a novelty in the current knowledge (S1 Table ).In fact, although their presence has been already described for different somatic tissues [42] and related to epigenetic reprogramming and gene expression regulation [3,43], piRNAs specifically located in mitochondria are only reported in mammalian cancer cells [22,44] and inferred by bioinformatic analyses on RNA extracted from the whole cell lysates derived from murine primordial germ cells, spermatozoa oocytes, zygotes and embryonic gonadal cells [21,41].Two studies have demonstrated a role of the mitochondrial protein MitoPLD/Zucchini (mammalian/Drosophila homologues) in piRNA biogenesis and piRNA mediated silencing in mouse and fly germlines [45,46]; finally, PIWIL-1, a human homologue of the mouse MIWI protein, has been detected in the mitochondria of human cells [44] pointing for a mitochondrial role in the biogenesis of piRNAs.Our discovery of piRNA populations in mitochondria of somatic tissues strongly supports the idea that the organelles may be an important and general player in piRNA pathways, beyond germ line tissues.Intriguingly, their expression pattern appears different for analyzed tissues, showing their massive presence in testis and some overlapping among analyzed tissues (Fig 3).Analogously to the so-called mitomiR, those piRNAs can be indicated as mitopiRNAs.

mitomiRNAs
The RNAseq analyses revealed similar percentages of miRNAs on the annotated ones in the mitochondria of gastrocnemius muscle, liver and WAT and a higher content in testis (Fig 2).
The overall expression pattern of mitomiRs is distinctive for the different tissues and ranges from few read counts to more than 1000, whose frequency counts is smaller compared to others (Fig 4A and 4B).We hypothesized that some miRNAs found in the mitochondria could be related to the basic functions of the organelle, whereas others could be mainly associated to specialized functions depending on the tissue.To explore that hypothesis, all the miRNAs found in the mitochondria of each tissue, from very low (> 1) to higher read counts were analyzed, resulting in a core of 126 shared mitomiRs, some others overlapping among the tissues and a number of tissue-specific mitomiRs (Fig 4C and S2 Table).Different papers and related datasets consistently report their expression also in tissues [47][48][49], and for some of them targets related to mitochondria functioning have been reported (Tables 1-5), thus further connecting their relevance for mitochondria biology.
It would be interesting to address the question whether miRNAs found in mitochondria reflect the general content of the cells or are differentially distributed between cytosol and mitochondria.To gain some insights into the issue, we quantify the level of specific miRNAs in cytosol and mitochondria fractions; in particular, we analyzed the level of the 1 st and the 10 th of tissues-shared top-10 mitomiRs according to their reads count mean, i.e. miR-122-5p and miR-133a-3p (Table 1), and two other miRNAs, miR-29a and miR-15b, respectively 12and 17-fold less expressed in comparison to miR-133a (S2 Table ).We found that in most cases microRNAs are differentially distributed between cytosol and mitochondria, appearing significantly enriched in mitochondria, (e.g., miR-122 in gastrocnemius, miR-29a in liver, miR-122 and miR-15b in testis, miR-133a in WAT), or in lower amount in comparison to the cytosol (e.g., miR-133a, miR-29a and miR-15b in gastrocnemius, miR-122 in liver, miR-29a in testis and WAT, miR-15b in WAT), pointing to a potential relevance of selective translocation into the organelles (Fig 5).

mitomiR targetomes and pathways
miRNAs exert their function by down-regulating multiple target transcripts, since their pairing can tolerate mismatches, beyond the perfect complementarity required for miRNA seed sequence (nucleotides 2-8).By screening the published papers on PubMed (searching keywords: "miRNA name AND mitochondria"), we found that various mitomiRs detected in our analyses have experimentally validated targets related to mitochondria function (Tables 1-5); however, gaps in the knowledge emerge for many others.To gain insights into biological functions of identified mitomiRs, the whole targetomes of all tissue-shared and all tissue-specific mitomiRs and related biological pathways were analyzed.In brief, we obtained a list of 4840 non-redundant transcripts experimentally validated or predicted by both TargetScan and miRDB tools as targets of the 126 tissue-shared mitomiRs; a list of 705 non redundant transcripts for the 21 mitomiR specifically found in gastrocnemius muscle; a list of 632 for the 30 mitomiRs for liver; a list of 4897 for the 174 testis-specific mitomiRs; a list of 520 for 18 WATspecific mitomiRs.
Enrichment of biological pathways supplied by KEGG on the targetome of tissue-shared mitomiRs shows an overrepresentation of different signaling pathways among the top p-values scored, beyond "Pathways in cancer", the most studied field in relation to microRNAs (Fig 6 and S3 Table).That observation is also applicable on KEGG analysis performed on targetome of testis specific mitomiRs, although including some different signaling pathways.The enrichment of biological pathways implied in different signaling routes was probably expected for miRNAs associated to mitochondria, since the central role of the organelle in metabolism and homeostasis and the results here reported highlight the potential strong miRNA contribution.Behind those pathways enriched for mitomiRs arising from all analyzed tissues, others are specifically found for tissue-specific mitomiRs; e.g, "insuline resistance" and "Type II diabetes mellitus" pathways were found for liver-specific mitomiRs; "Wnt signling pathway", wellknown to be related to gluconeogenesis, was found for WAT; "Regulation of actin cytoskeleton" pathway emerges significantly enriched for mitomiRs in testis, where actin remodeling contributes to proper spermatogenesis in terms of morphology and then motility [74].Overall, those results support the idea that regulation via miRNAs could be functionally relevant for mitochondria biology, since a pool of miRNAs could be related to the basic function of mitochondria shared by different cell types, and others may be related to the specialized functions of each tissue.
KEGG annotations was helpful to put mitochondrial pathways potentially governed by miRNAs in a larger context, but we were also interested in exploring mitochondria-centric pathway annotations complementary to the broader pathway database.To that aim, the targetome lists were matched with 1140 transcripts encoding proteins with mitochondrial localization in the database of MitoCarta3.0, manually curated by adding experimentally validated transcripts targets as resulting from PubMed for our detected mitomiRs (Fig 7 and S4 Table ).The assignment to the "MitoPathways" spanning the 7 categories related to mitochondria Table 1.Tissues-shared top-10 mitomiRs according to their reads count.

Discussion
Mitochondria are "multitasked" organelles, tailoring even their morphology, number and distribution to execute specialized functions in different cell types and/or different environments, behind the oxidative phosphorylation.Proteins and protein-driven pathways involved in mitochondria homeostasis are well known, whereas non-coding RNA contribution is still overlooked.In order to gain insights into possible new and specialized function of small RNAs related to mitochondria biology, here we performed a survey of small RNA content of mitochondria purified from different tissues that particularly relies on mitochondria functioning (Fig 1).
It is described for the first time the presence of mitopiRNAs in somatic tissues, behind testis tissues, thus supporting the idea of a strong and ubiquitous connection between mitochondria biology and piRNA pathways (Figs 2 and 3).piRNAs represents an emerging and fascinating class of short RNA, longer than miRNAs, usually 24-35 nt, whose biogenesis and functioning still require further investigation along with validation of those piRNAs identified only on their size, rather than on their PIWI-interaction or functional roles.It is well established that piRNAs are enriched in germ cells, where they modulate repression of transposable elements, and probably also gene expression regulation [5][6][7].With regard to biogenesis, piRNAs are generally transcribed as longer precursors from specific genomic regions termed "piRNA clusters"; then, they are processed into mature piRNAs trough two interconnected pathways, ping- pong amplification mediated by PIWI proteins in ribonucleoprotein granules, and phasing, where PIWI-loaded piRNA precursors are translocated to the mitochondrial outer membrane by RNA helicases and further cleaved by endonuclease Zucchini/PLD [75].Also PIWIL-1 was observed to localize to mitochondria, thus raising the possibility of a mitochondrial role in the biogenesis of piRNAs [44].piRNAs in mitochondria have been already reported, but only for mammalian cancer cells [22,44] and inferred by bioinformatic analyses on RNA purified from the whole cell lysates deriving from murine primordial germ cells, spermatozoa, oocytes, zygotes and embryonic gonadal cells [21,41].Those studies also suggest the possibility that some piRNAs may be generated also from mitochondrial genome.In the mitochondria purified from the tissues, we found different populations of piRNAs, shared or tissue-specific, but none of them seems to be produced by mitochondrial genome.To assess that, all piRNA sequences contained in piRBAse were aligned to mouse mitochondrial genome, finding 7.29% of perfectly matching piRNA full sequences (S5 Table ), and thus possibly deriving from mitochondrial genome, consistently with data already reported [41]; however, no reads related to those perfectly matching piRNAs were found in our samples.Different reads related to the other piRNAs, but not completely covering their full sequences, cross-map to mitochondrial genome (S6 Table and S1 Fig), probably due to the presence of NUMTs, nuclear-mitochondrial identical DNA sequences, found in many metazoans [14,76].Overall, our data indicate that mitopiRNAs populations found here and commonly or specifically expressed in the different tissues are generated by nuclear genome, translocated to mitochondria and involved in a piRNA-mediated communication between mitochondria and nucleus, whose mechanisms and biological meaning deserves further investigation.
With regard to miRNAs, it is known that they can regulate different pathways related to mitochondria, e.g. by targeting transcripts encoding proteins involved in TCA cycle, OXPHOS, and mitochondrial dynamics [13,75].Although miRNAs are largely operative in the cytoplasm, some other papers report their localization in the nucleus and mitochondria [13,[15][16][17][18][19][20][21][22], whose significance and mechanism of action are still largely unknown, since the debated presence of RISC components inside the organelle.Here we describe a mitomiRs landscape, finding a group of mitomiRs shared by all tissues analyzed, and thus potentially endowed with a regulatory role required for general mitochondria functioning, and tissue-specific mitomiRs, potentially related to specialized functions dependent on cell context (Fig 4).Among the top-10 expressed miRNAs, different ones have experimentally validated targets related to mitochondria functioning, as reported in Tables 1-5; however, for many of them studies demonstrating possible functional links with mitochondria are still missing and could be inspired by their finding inside the organelle.To gain insights into biological functions of detected mitomiRs, we extended our analyses to the whole targetomes of all tissue-shared and all tissue-specific mitomiRs, considering both experimentally validated and predicted targets and their enrichment in biological pathways.KEGG analyses highlighted that the tissue-shared mitomiRs greatly contribute to the central role of mitochondria in different signaling routes, overrepresented among the top p-values scored pathways; mitomiRs detected in specific tissues potentially contribute to regulate biological pathways typical of those tissues, as detailed in Results, overall supporting the idea of a shared population of mitomiRs for basic functioning and others for specialized functions related to the different tissues.In a more mitochondria-centric view, exploitation of MitoCarta3.0 database further confirm the presence of targets encoding proteins related to mitochondria functioning in the whole targetomes and their involvement in pathways specific of the organelle such as "Metabolism" and "OXPHOS" ( Fig 7 and S4 Table).
With regard to mapping of mitomiRs to nuclear genome, intriguingly, 8 out of 10 most expressed mitomiRs specifically found in testis is produced by X chromosome (Table 5), consistently with data already reported and arguing a special role of miRNAs from X chromosome in male fertility [77,78].
Finally, it is debated if mitomiRs could also derive from mitochondrial genome, although the mainstream view is that they originate from the nucleus and then translocated to mitochondria as illlustrated in Introduction.In order to give our contribution to the debate, we followed the same workflow described above for piRNAs to assess the miRNAs origin in our sample: all murine miRNAs deposited in miRbase were aligned to mitochondrial genome, finding only mmu-miR-6481 and mmu-miR-12206-5p sequences perfectly matching as full sequences to mitochondrial genome.However, those 2 miRNAs are not represented in our RNAseq data, instead finding reads related to other miRNAs, although not covering their full sequences (S6 Table and S1 Fig).Also taking into consideration the absence of Dicer in mitochondria [18,47,79], the data suggest that in our samples the reads related to miRNAs and mapping to mitochondrial genome can be attributed to miRNAs produced by their wellknown chromosomal location, and could cross-map to mitochondrial genome for the presence of NUMTs, as also discussed above for mitopiRNAs [14,76].Besides miRNAs and piR-NAs, other small non-coding RNA biotypes (generally called mitosRNA) could be produced by mitochondrial genome, not affected by Dicer absence, but probably products of currently unidentified mitochondrial ribonucleases, and whose biogenesis and role still remain unknown requiring further investigation [47].
Overall, our analyses represent a large data platform obtained by RNAseq on highly purified mitochondria from different tissues, suggesting the idea of a "basic" small RNA population residing into mitochondria and shared by different tissues, and others restricted to specific tissues.The significantly different profiles between mitochondria and cytosol fractions for some analyzed mitomiRs (Fig 5) support the relevance of miRNAs for mitochondria functioning and a potential selective translocation.The next challenge will be the planning of functional studies to deepen our understanding of small RNA processing, the shuttling among different subcellular compartments and their governed regulatory networks contributing to mitochondria biology.
Fig 1).The percentage of detected mitochondrial smallRNA populations on the annotated ones for each tissue are reported in Fig 2 and detailed below.
mitopiRNAs piRNAs are largely present in the mitochondria of testis, representing the most predominant population of small RNAs (Fig 2)

Fig 3 .
Fig 3. mitopiRNAs expression analyses.(A) Normalized expression unsupervised heatmap relative to piRNAs in the mitochondria of the indicated tissues, obtained by using the software ComplexHeatmap, a package of R. (B) The Venn diagram displaying the numbers of overlapping mitopiRNAs (Raw Count Mean �1) among the different tissues.https://doi.org/10.1371/journal.pone.0293644.g003

Fig 4 .
Fig 4. mitomiRNAs expression analyses.(A) Normalized expression unsupervised heatmaps relative to mitomiRs of the indicated tissues, obtained by using the software ComplexHeatmap, a package of R. (B) Distribution of mitomiR levels with respect to frequency; numbers of sequence reads (showed on the bars) were taken as miRNA levels and the values are represented in the form of range of values in mitochondrial sRNA libraries of each tissue.(C) The Venn diagram displaying the numbers of overlapping mitomiRs among the different tissues.https://doi.org/10.1371/journal.pone.0293644.g004

Fig 5 .Fig 6 .
Fig 5. Expression patterns of miRNAs in cytosol and mitochondria fractions from the different tissues.Expression analyses of the indicated miRNAs were performed by using Taqman assays on the RNA extracted from cytosol and mitochondria fractions of the indicated tissues.All analyses were performed at least in triplicate on each sample and mean + sd is reported.p-values at Student's t-test were *p < 0.05, **p < 0.01, or ***p < 0.001.https://doi.org/10.1371/journal.pone.0293644.g005